# # Pre-analysis plan for survey experiment 
#  	implemented in German Internet Panel 
#  	wave 24 (July 2016). 
# 
#
# # This version: 2016/07/11

setwd("")
set.seed(42)

library(car)
library(classInt)
library(dplyr)
library(ggplot2)
library(lmtest)
library(readstata13)
library(rstanarm)
library(sandwich)
library(stringr)
library(texreg)
library(tidyr)

get_ols <- function(f,d){
	m <- lm(f, data=d)
	m <- coeftest(m, vcov = vcovHC(m, "HC3"))
	return(m)
	}


# Recode Strings 
age_recode <- 
"'1. 16-19'='16-19';
'2. 20-24'='20-24';
'3. 25-29'='25-29';
'4. 30-34'='30-34';
'5. 35-39'='35-39';
'6. 40-44'='40-44';
'7. 45-49'='45-49';
'8. 50-54'='50-54';
'9. 55-59'='55-59';
'10. 60-64'='60-64';
'11. 65-69'='65-69';
'12. 70-74'='70-74';
'13. 75-79'='75-79';
else=NA"

gender_recode <- 
"'1. mÃ¤nnlich'='male'; 
'2. weiblich'='female';
else=NA"

number_hh <- 
"'1. 1'='1';
'2. 2'='2';
'3. 3'='3';
'4. 4'='4';
'5. 5'='5';
'6. 6 und mehr'='6';
else=NA"

lr_recode <- "'1. links'='1';
'2'='2';
'3'='3';
'4'='4';
'5'='5';
'6'='6';
'7'='7';
'8'='8';
'9'='9';
'10'='10';
'11. rechts'='11';
'-99. weiß nicht/keine Meinung'='0';
else=NA"

occup_recode <- 
"'1. Vollzeiterwerbstätig'='1';
'2. Teilzeiterwerbstätig'='1';
'3. Altersteilzeit (unabhängig davon, ob in der Arbeits- oder Freistellungsphase bef'='1';
'4. Geringfügig erwerbstätig, 400-Euro-Job, Minijob'='1';
'5. „Ein-Euro-Job“ (bei Bezug von Arbeitslosengeld II)'='0';
'6. Gelegentlich oder unregelmäßig beschäftigt'='0';
'7. In einer beruflichen Ausbildung/Lehre'='0';
'8. In Umschulung'='0';
'9. Freiwilliger Wehrdienst, Bundesfreiwilligendienst'='0';
'10. Freiwilliges Soziales/Ökologisches/Kulturelles Jahr'='0';
'11. Mutterschafts-, Erziehungsurlaub, Elternzeit oder sonstige Beurlaubung (Alterste'='1';
'12. Schüler/-in an einer allgemeinbildenden Schule'='0';
'13. Student/-in'='0';
'14. Rentner/-in, Pensionär/-in, im Vorruhestand'='0';
'15. Arbeitslos'='0';
'16. Dauerhaft erwerbsunfähig'='0';
'17. Hausfrau/Hausmann'='0';
else=NA"



################
################	
################
# 	MAKE DATA     #
################
################
################


# Demographics from wave 19
######################

w19 <- read.dta13("ZA5985_v1-0-0.dta", 
	nonint.factors=TRUE, generate.factors=TRUE, encoding="UTF-8")

w19 <- w19 %>% select(hhid_g, gender_15, age_cat_15, educ_school_15, educ_job_15,
	number_hh_members_15, occupation_15, state, AA19040a)

w19 <- w19 %>% 
	mutate( gender_15 = recode(gender_15, gender_recode),
	 	   age_cat_15 = recode(age_cat_15, age_recode),
	 	   number_hh_members_15 = recode(number_hh_members_15, number_hh, 
	 	   	as.factor.result=FALSE),
	 	   AA19040a = recode(AA19040a, lr_recode, as.factor.result=FALSE),
	 	   state=unclass(state)-6,
	 	   edu = 1, educ_school_15 = unclass(educ_school_15),
	 	   educ_job_15 = unclass(educ_job_15),
	 	   edu = ifelse(educ_school_15 %in% c(11,10,12), 2, edu),
	 	   edu = ifelse(educ_job_15 %in% c(17,16,15,14), 3, edu),
	 	   occupation_15 = recode(occupation_15, occup_recode) ) %>% 
	rename( lr = AA19040a) %>% 
	select(-educ_job_15,-educ_school_15)



# Exp Conditions (fake)
#################

N <- nrow(w19)

exp_rule <- sample(c(1,2,3),N, replace=TRUE)
w19$exp_rule <- factor(exp_rule, labels=c("maj", "super", "unan"))

exp_decision <- sample(c(1,2),N, replace=TRUE)
w19$exp_decision <- factor(exp_decision, labels=c("adop", "rej"))

exp_ideology  <- sample(c(1,2),N, replace=TRUE)
w19$exp_ideology <- factor(exp_ideology, labels=c("left", "right"))

exp_gender <- sample(c(1,2),N, replace=TRUE)
w19$exp_gender <- factor(exp_gender, labels=c("male", "female"))

exp_experience <- sample(c(1,2,3),N, replace=TRUE)
w19$exp_experience <- factor(exp_experience, labels=c("fresh", "jun", "senior"))

exp_membership <- sample(c(1,2),N, replace=TRUE)
w19$exp_membership <- factor(exp_membership, labels=c("edu", "econ"))


# Exp Outcomes (fake)
################

w19$exp_prr_blf_vtg <- runif(N, 0, 100)
w19$exp_ptrr_blf_vtg <- runif(N, 0, 100)

w19$outcome <- w19$exp_prr_blf_vtg-w19$exp_ptrr_blf_vtg

w19 <- w19 %>% 
	mutate( exp_prr_blf_vtg_group=cut(exp_prr_blf_vtg, 
			c(0,25,50,75,100), labels=FALSE) )

response_time <- runif(N, 0, 60)
response_time_brks <- with(classIntervals(response_time, n=3, style="quantile"), brks)
w19$time <- cut(response_time, response_time_brks, labels=FALSE)

gip <- w19



################
################	
################
# Theory-Data Plot #
###############
################
################


Rs <- c(8,11,14)
M <- 14
qs <- seq(0.10,0.90, length=100)

I <- 3
J <- length(qs)

beliefA <- beliefB <- matrix(NA, I, J)

for(i in 1:I){
	for(j in 1:J){
	
	x = 1-pbinom(Rs[i]-2,M-1,qs[j]) 
	y = 1-pbinom(Rs[i]-1,M-1,qs[j]) 

	# Adoption 
	numA <- qs[j]*x
	denA <- qs[j]*x + (1-qs[j])*y
	beliefA[i,j] <- numA/denA 

	# Rejection
	numB <- qs[j]*(1-x)
	denB <- numB + (1-qs[j])*(1-y)  
	beliefB[i,j] <- numB/denB 
	}
}

beliefA <- beliefA*100
beliefB <- beliefB*100
qs <- qs * 100

theory <- data.frame(qs, t(beliefA), t(beliefB) )
colnames(theory) <- c("prr_blf", "ptrr_blf|maj,adop", 
	"ptrr_blf|super,adop", "ptrr_blf|unan,adop",
	"ptrr_blf|maj,rej", "ptrr_blf|super,rej", 
	"ptrr_blf|unan,rej")

theory <- theory %>% gather(cond, ptrr_blf, -prr_blf) %>% 
		separate(cond, c("exp_rule", "exp_decision"), ",") %>% 
		mutate(exp_rule = str_replace(exp_rule, "ptrr_blf\\|", ""))

ggplot(theory, aes(prr_blf, ptrr_blf,color=exp_rule)) + 
	facet_grid(. ~ exp_decision) + 
	geom_line(lwd=1.5) + geom_point(aes(exp_prr_blf_vtg, 
		exp_ptrr_blf_vtg), data=gip) +
	ylab("Posterior Probability") + 
	xlab("Prior Probability")  



###########
###########
###########
# ESTIMATES #
###########
###########
###########

gip <- na.omit(gip)

# Slice Data 
adop <- filter(gip, exp_decision=="adop")
rej <- filter(gip, exp_decision=="rej")


######
# ols  #
######

f <- outcome ~ factor(exp_rule)

f_full <- update(f, . ~ . + factor(edu) + 
	factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

f_full_edu <- update(f, . ~ .  * factor(edu) + 
	factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

f_full_time <- update(f, . ~ .  * factor(time) + 
	factor(edu) + factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

f_full_prior <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) + 
	factor(edu) + factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

f_full_prior_edu <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) * factor(edu) + 
	factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

f_full_prior_time <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) * factor(time) + 
	factor(edu) + 
	factor(gender_15) + factor(age_cat_15) + 
	factor(occupation_15) + factor(lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) )

formula_batch <- list(f,
	f_full,
	f_full_edu,
	f_full_time,
	f_full_prior,
	f_full_prior_edu,
	f_full_prior_time)

ols_adop <- lapply(formula_batch, function(formula) get_ols(formula, adop) )
ols_rej <- lapply(formula_batch, function(formula) get_ols(formula, rej) )




###########
# stan_lmer  #
###########

f <- outcome ~ factor(exp_rule) + (1|hhid_g) 

f_full <- update(f, . ~ . + factor(edu) + 
	factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) + 
	(1|hhid_g)  )

f_full_edu <- update(f, . ~ .  * factor(edu) + 
	factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) + 
	(1|hhid_g)  )

f_full_time <- update(f, . ~ .  * factor(time) + 
	factor(edu) + factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) + 
	(1|hhid_g)  )

f_full_prior <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) + 
	factor(edu) + factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) +
	(1|hhid_g)  )

f_full_prior_edu <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) * factor(edu) + 
	factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) + 
	(1|hhid_g) )

f_full_prior_time <- update(f, . ~ .  * factor(exp_prr_blf_vtg_group) * factor(time) + 
	factor(edu) + 
	factor(gender_15) + (1|age_cat_15) + 
	factor(occupation_15) + (1|lr) + 
	factor(exp_ideology) +
	factor(exp_gender) + 
	factor(exp_experience) +
	factor(exp_membership) + 
	(1|hhid_g) )

formula_batch <- list(f,
	f_full,
	f_full_edu,
	f_full_time,
	f_full_prior,
	f_full_prior_edu,
	f_full_prior_time)


prior <- normal(0,2.5)
prior_intercept <- normal(0,50)

lmer_adop <- lapply(formula_batch, function(formula) stan_lmer(formula, 
	data=adop, prior=prior, prior_intercept=prior_intercept) )
lmer_rej <- lapply(formula_batch, function(formula) stan_lmer(formula, 
	data=rej, prior=prior, prior_intercept=prior_intercept) )


#######
# more  #
#######

# - we estimate all models also without Speeders which are defined as all respondents who 
#	finished the entire survey in less than half of the median time.


